% COZZI AND IMPULLITTI (2015), GLOBALIZATION AND WAGE POLARIZATION, REVIEW
% OF ECONOMICS AND STATISTICS

% RUN FILE FOR COMPUTATION OF EQUILIBRIUM AND COMPARATIVE STATICS

clear all

format short

x0plus=[ 0.84542445670327   0.84266619012432   1];

global sig1 gam beta phi AD tau epsi fi1 n lambda eta rho eserv fi zd1 zf1

AF0=  1:0.01:7.5; % Technology parameter innovation F. An increase in AF decreases the technology gap

%======================== Externally calibrated ===========================


rho=0.03;      % discount factor

lambda=1.3;    % quality jump
n=0.0114;      % population growth rate
gam = 0.6;     % upper bound for share of skilled workers
eta = 1.5;     % Elasticity of substitution
fi1 = 0.0256;  % R&D difficulty parameter
V = 52;        % Life duration after college
TH = 4;        % Year of college
tau =1.74;     % Trade cost
zf1=1;         % Technology parameter production F

% ================  Internally calibrated  ================================

par_vec = [ 1.512786682942995   0.433971610241538   6.582402166785685   0.512486624266001 0.023413083220775   5.089327896809821];

zd1= par_vec(1);       % Technology parameter production D
eserv = par_vec(2);    % Service efficiency parameter
epsi = par_vec(3);     % Ability distribution curvature
beta= par_vec(4);      % Factor intensity production
phi = par_vec(5);      % Factor intensity innovation
AD = par_vec(6);       % Technology parameter innovation D
   

ntau=length(AF0);


sig1=((exp(rho*V)-1)/(exp(rho*(V-TH))-1));
fi = ((exp(n*(V-TH))-1)/(exp(n*V)-1));


options=optimset('MaxFunEvals',100000000,'MaxIter',100000000,'Tolfun',0.0000000001,'TolX',0.00000000001);

XX01=zeros(ntau,1);
XX02=zeros(ntau,1);
XX03=zeros(ntau,1);
XX1=zeros(ntau,1);
XX2=zeros(ntau,1);
XX3=zeros(ntau,1);
XX4=zeros(ntau,1);
XX5=zeros(ntau,1); 
XX6=zeros(ntau,1);
XX7=zeros(ntau,1);
XX8=zeros(ntau,1);
XX9=zeros(ntau,1);
XX10=zeros(ntau,1);
XX11=zeros(ntau,1);
XX12=zeros(ntau,1);
XX13=zeros(ntau,1);
XX14=zeros(ntau,1);
XX15=zeros(ntau,1);
XX16=zeros(ntau,1);
XX16b=zeros(ntau,1);
XX17=zeros(ntau,1);
XX17b=zeros(ntau,1);
XX17plot=zeros(ntau,1);
XX17bplot=zeros(ntau,1);
XX18=zeros(ntau,1);
XX19=zeros(ntau,1);
XX19b=zeros(ntau,1);
XX20=zeros(ntau,1);
XX21=zeros(ntau,1);
XX22=zeros(ntau,1);
XX23=zeros(ntau,1);
XX24=zeros(ntau,1);
XX25=zeros(ntau,1);
XX26=zeros(ntau,1);
XX27=zeros(ntau,1);
XX28=zeros(ntau,1);
XX29=zeros(ntau,1);
XX30=zeros(ntau,1);
XX31=zeros(ntau,1);
XX32=zeros(ntau,1);


for i=1:ntau
    AF=AF0(i);

 [x,fval,exitflag,output]=fsolve(@CI_compstat,x0plus,options,AF);

 
 exitflag

            if exitflag<0
                disp('Problem here!');
                pause
            return
            end
            

the0D = x(1);
the0F = x(2);
wDL = x(3);

wDH = (sig1*(the0D/(the0D-gam))*x(3));
wFL = (((the0F-gam)/(sig1*the0F))^(1-beta));
wFH = (sig1*(the0F/(the0F-gam))*wFL);

% ======================== TECHNOLOGY  ====================================

ZD = (1/zd1)*(wDL^beta)*(wDH^(1-beta));
ZF = (1/zf1)*(wFL^beta)*(wFH^(1-beta));
FD = (1/AD)*((wDL^phi)*(wDH)^(1-phi));
FF = (1/AF)*((wFL^phi)*(wFH)^(1-phi));

   % Marginal products: production
   
ZDL = (1/zd1)*(beta*((wDH/wDL)^(1-beta)));
ZDH = (1/zd1)*((1-beta)*((wDH/wDL)^(-beta)));
ZFL = (1/zf1)*(beta*((wFH/wFL)^(1-beta)));
ZFH = (1/zf1)*((1-beta)*((wFH/wFL)^(-beta)));
   
   % Marginal products: R&D
   
FDL = (phi*(1/AD)*((wDH/wDL)^(1-phi)));
FDH = ((1-phi)*(1/AD)*((wDH/wDL)^(-phi)));
FFL = (phi*(1/AF)*((wFH/wFL)^(1-phi)));
FFH = ((1-phi)*(1/AF)*((wFH/wFL)^(-phi)));

% =========================== CUTOFFS =====================================

   % Skill cutoff to hire service workers
   
theHSD_fun = @(y) 1-y^epsi-(((wDH/wDL)*(1-eserv)*(y-gam))^epsi);
theHSF_fun = @(y) 1-y^epsi-(((wFH/wFL)*(1-eserv)*(y-gam))^epsi);

theHSD0 = 0.6;
theHSF0 = 0.6;

theHSD = csolve(theHSD_fun,theHSD0,[],1e-7,300);
theHSF = csolve(theHSF_fun,theHSF0,[],1e-7,300);

   % Skill cutoff for unskilled in production

theLSD = ((1-(theHSD^epsi))^(1/epsi));
theLSF = ((1-(theHSF^epsi))^(1/epsi));
  
   % Labor supplies
   
LDm = ((epsi/(1+epsi))*(the0D^(1+epsi)-theLSD^(1+epsi)));  % unskilled supply of production workers
LFm = ((epsi/(1+epsi))*(the0F^(1+epsi)-theLSF^(1+epsi)));
HD = ((1-the0D^epsi)*fi*(((epsi*((1-theHSD^(1+epsi))*(2 - eserv) + theHSD^(1+epsi) - the0D^(1+epsi)))/((1+epsi)*(1-the0D^epsi)))-gam)); % Skilled labor supply
HF = ((1-the0F^epsi)*fi*(((epsi*((1-theHSF^(1+epsi))*(2 - eserv) + theHSF^(1+epsi) - the0F^(1+epsi)))/((1+epsi)*(1-the0F^epsi)))-gam));

   % Average workers skills
   
theDm_avg = ((epsi/(1+epsi)) *((the0D)^(1+epsi)-(theLSD)^(1+epsi))/(the0D^epsi-theLSD^epsi)); % Unkilled: Production
theFm_avg = ((epsi/(1+epsi)) *((the0F)^(1+epsi)-(theLSF)^(1+epsi))/(the0F^epsi-theLSF^epsi));
theDS_avg = ((epsi/(1+epsi)) *theLSD);  % Unskilled: Services corrected
theFS_avg = ((epsi/(1+epsi)) *theLSF);
theDH_avg = (((epsi*((1-theHSD^(1+epsi)) + theHSD^(1+epsi) - the0D^(1+epsi)))/((1+epsi)*(1-the0D^epsi)))-gam);  % This average corrects for the increased supply of labour by skills: it is the average over the individuals of the averages per unit time, which means that those who work 2-eserv hours have their wages divided by (2-eserv). This is why we do not see (2-eserv) there, which of course we see it in the aggregate skill supply
theFH_avg = (((epsi*((1-theHSF^(1+epsi)) + theHSF^(1+epsi) - the0F^(1+epsi)))/((1+epsi)*(1-the0F^epsi)))-gam);  % This average corrects for the increased supply of labour by skills: it is the average over the individuals of the averages per unit time, which means that those who work 2-eserv hours have their wages divided by (2-eserv). This is why we do not see (2-eserv) there, which of course we see it in the aggregate skill supply
theDunskill_avg = ((epsi/(1+epsi))*the0D); % Average unskilled (service+production)
theFunskill_avg = ((epsi/(1+epsi))*the0F);

% ======================= OTHER ENDOGENOUS VARIABLES ======================

% Here we report some expressions we use the reduced the original
% equilibrium system to 3 equations.

qF_til = ((1-fi1)*((lambda^(eta-1))-1)/n);
Fr = (FD/FF);
cD_til =( (Fr*(lambda-1)-(lambda-ZD*tau)) / ((lambda-ZD/tau)-(lambda-1)*Fr) );
x_til = ((1/FF)* (((lambda-1)/lambda)*(1+cD_til))/ (rho+(n/((1-fi1)*((lambda^(eta-1))-1))) + (n*fi1/(1-fi1))));
IF = (1/qF_til)*(1/(cD_til+1));
CF = ((1/IF)*(LFm / (qF_til*(ZFL/lambda)*(1+cD_til) + x_til*FFL)));
qF = (qF_til*IF);
qD = (1-(qF));
ID = ((1/qF_til)-IF);
CD = (qD/qF)*CF;
xindex = ((x_til)*CF); % R&D difficulty indiex
CapitalD = FD*xindex*qD;% total value of the domestic firms
CapitalF = FF*xindex*qF;% total value of the domestic firms
GDPD = qD*(CD+CF) + (wDL*FDL + wDH*FDH)*xindex*ID; % THE WAGES OF THE SERVICE WORKERS ARE PAID BY THE SKILLED WORKERS, SO THEY CANCEL OUT FROM AGGREGATE INCOME!  
GDPF = qF*(CD+CF) + (FFL + wFH*FFH)*xindex*IF; %THE WAGES OF THE SERVICE WORKERS ARE PAID BY THE SKILLED WORKERS, SO THEY CANCEL OUT FROM AGGREGATE INCOME! 
betaPikD = CapitalD/GDPD; % Wealth-income ratio for country D
betaPikF = CapitalF/GDPF; % Wealth-income ratio for country F

% ======================== Output variables ===============================

XX00(i,1) = FD;
XX0F(i,1)= FF;
XXZF(i,1)= ZF;
XX01(i,1) = the0D;
XX02(i,1) = the0F;
XX03(i,1) = wDL;
XX1(i,1) = CD;
XX2(i,1) = CF;
XX3(i,1) = ID;
XX4 (i,1) = IF;
XX5 (i,1) = wDL;
XX6 (i,1) = wDH;
XX7 (i,1) = wFL;
XX8 (i,1) = wFH;
XX9 (i,1) = ZD;
XX10 (i,1) = wDH*theDH_avg; %avg skilled wages D
XX11 (i,1) = wFH*theFH_avg; %avg skilled wages F
XX12 (i,1) = wDL*theDm_avg;  %avg production wages D
XX13 (i,1) = wFL*theFm_avg;  %avg production wages F
XX14 (i,1) = wDL*theDS_avg; % avg service wages D
XX15 (i,1) = wFL*theFS_avg; % avg service wages F
XX16 (i,1) = (wDH*theDH_avg)/(wDL*theDm_avg); % Avg. Skilled / Avg. Unskilled Wage D
XX16b (i,1) = (wDH*theDH_avg)/(wDL*theDunskill_avg ); % skill premium D
XX17 (i,1) = (wFH*theFH_avg)/(wFL*theFm_avg); % Avg. Skilled / Avg. Unskilled Wage F
XX17b (i,1) = (wFH*theFH_avg)/(wFL*theFunskill_avg ); % skill premium F
XX17plot(i,1) = 1.462;
XX17bplot(i,1) = 1.51;
XX18 (i,1) = (theDm_avg)/(theDS_avg); % Production premium D
XX19 (i,1) = (theFm_avg)/(theFS_avg);  % Production premium F
XX19b(i,1) = 1.35;
XX20 (i,1) = (1-((the0D)^epsi))*fi;  % share of skilled workers (physical units, not accounting for efficiency)
XX21 (i,1) = ((the0D)^epsi - theLSD^epsi); % share of unskilled (no service)
XX22 (i,1) = ((theLSD)^epsi); % share of workers in service
GDP = (qD*((CD/ZD + CF*tau)*(ZD/lambda) + (rho-n)*FD)); 
XX24 (i,1) = ID/(ID+IF); % Patent share
XX25 (i,1) = xindex;
XX26(i,1) = (1-the0D^epsi)*((exp(n*V)-exp(n*(V-TH)))/(exp(n*V)-1)); %share of students D
XX27(i,1) = CapitalD;
XX28(i,1) = CapitalF;
XX29(i,1) = GDPD;
XX30(i,1) = GDPF;
XX31(i,1) = betaPikD;
XX32(i,1) = betaPikF;

x0plus=x;

end


disp('Initial equilibrium values: the0D the0F WDL') 

init = [XX01(1,1) XX02(1,1) XX03(1,1)]

%========================== Comparative static results  ===================

disp('INITIAL MOMENTS')
disp(' Skill/Prod_wage    Skill Prem.       Prod/Service_wage       % Skilled        % Production   % Service   Patent share  Beta Piketty')
MomentsInitial =  [XX16(1,1), XX16b(1,1), XX18(1,1), XX20(1,1), XX21(1,1), XX22(1,1), XX24(1,1), XX31(1,1)]

disp('FINAL MOMENTS')
disp(' Skill/Prod_wage    Skill Prem.       Prod/Service_wage       % Skilled        % Production   % Service   Patent share  Beta Piketty')
MomentsFinal =  [XX16(ntau,1), XX16b(ntau,1), XX18(ntau,1), XX20(ntau,1), XX21(ntau,1), XX22(ntau,1), XX24(ntau,1), XX31(ntau,1)]


% Share of variation in the data explained by the model

DSkilProdD =((XX16(ntau,1) - XX16(1,1) )/XX16(1,1))
DSkillPD =((XX16b(ntau,1) - XX16b(1,1) )/XX16b(1,1))
DUnskillpD =((XX18(ntau,1) - XX18(1,1) )/XX18(1,1))
DshareSkill =((XX20(ntau,1) - XX20(1,1) )/XX20(1,1))
DshareProduction =((XX21(ntau,1) - XX21(1,1) )/XX21(1,1))
DshareService =((XX22(ntau,1) - XX22(1,1) )/XX22(1,1))
Dpatshare =((XX24(ntau,1) - XX24(1,1) )/XX24(1,1))
DbetaPikD = ((XX31(ntau,1) - XX31(1,1) )/XX31(1,1));

disp('Model Explanatory Power of Changes in data')

Shareschange_explain = [DSkilProdD/0.21, DSkillPD/0.174  DUnskillpD/0.16, DshareSkill/0.25, DshareProduction/0.166, DshareService/0.10  Dpatshare/0.156, DbetaPikD/((2.82-0.67)/0.67)]

% DbetaPikD in data is:((2.82-.67)/.67) % US Corporate Equity Wealth/Income was 67% in 1980 and 282% in 2000


% ========================= FIGURES PDF paper =============================

xSize = 19; 
ySize = 11*2;
thick = 1.8;
figure;

% suptitle('Figure 4. Technology Gap and Labor Market Polarization') 

subplot(4,3,1)
plot(AF0,XX16,'LineWidth',thick,'color','k');
title('Skilled/unskilled wage D','Interpreter','latex','FontSize',10)
axis tight
xlabel('Technology gap, ${A^{F}}/{A^{D}}$','Interpreter','latex','FontSize',10)
set(gca,'FontSize',8,'TickLabelInterpreter','latex')
set(gcf,'Units','centimeters','Position',[0 0 xSize ySize],'PaperUnits','centimeters' ...
     ,'PaperPosition',[0 0 xSize ySize],'PaperSize',[xSize-2.7 ySize-2],'PaperPositionMode','auto')

subplot(4,3,2)
plot(AF0,XX16b,'LineWidth',thick,'color','k');
title('Skill Premium D','Interpreter','latex','FontSize',10)
axis tight
xlabel('Technology gap, ${A^{F}}/{A^{D}}$','Interpreter','latex','FontSize',10)
set(gca,'FontSize',8,'TickLabelInterpreter','latex')
set(gcf,'Units','centimeters','Position',[0 0 xSize ySize],'PaperUnits','centimeters' ...
     ,'PaperPosition',[0 0 xSize ySize],'PaperSize',[xSize-2.7 ySize-2],'PaperPositionMode','auto')
 
subplot(4,3,3)
plot(AF0,XX18,'LineWidth',thick,'color','k');
title('Unskilled/Service Wage D','Interpreter','latex','FontSize',10)
axis tight
xlabel('Technology gap, ${A^{F}}/{A^{D}}$','Interpreter','latex','FontSize',10)
set(gca,'FontSize',8,'TickLabelInterpreter','latex')

 set(gcf,'Units','centimeters','Position',[0 0 xSize ySize],'PaperUnits','centimeters' ...
     ,'PaperPosition',[0 0 xSize ySize],'PaperSize',[xSize-2.7 ySize-2],'PaperPositionMode','auto')
 
subplot(4,3,4)
plot(AF0,XX24,'LineWidth',thick,'color','k');
title('Patent Share D','Interpreter','latex','FontSize',10)
axis tight
xlabel('Technology gap, ${A^{F}}/{A^{D}}$','Interpreter','latex','FontSize',10)
set(gca,'FontSize',8,'TickLabelInterpreter','latex')

 set(gcf,'Units','centimeters','Position',[0 0 xSize ySize],'PaperUnits','centimeters' ...
     ,'PaperPosition',[0 0 xSize ySize],'PaperSize',[xSize-2.7 ySize-2],'PaperPositionMode','auto')
 
 subplot(4,3,5)
plot(AF0,XX20,'LineWidth',thick,'color','k');
title('Skilled Workers D','Interpreter','latex','FontSize',10)
axis tight
xlabel('Technology gap, ${A^{F}}/{A^{D}}$','Interpreter','latex','FontSize',10)
set(gca,'FontSize',8,'TickLabelInterpreter','latex')

 set(gcf,'Units','centimeters','Position',[0 0 xSize ySize],'PaperUnits','centimeters' ...
     ,'PaperPosition',[0 0 xSize ySize],'PaperSize',[xSize-2.7 ySize-2],'PaperPositionMode','auto')
 
  subplot(4,3,6)
plot(AF0,XX21,'LineWidth',thick,'color','k');
title('Unskilled Workers D','Interpreter','latex','FontSize',10)
axis tight
xlabel('Technology gap, ${A^{F}}/{A^{D}}$','Interpreter','latex','FontSize',10)
set(gca,'FontSize',8,'TickLabelInterpreter','latex')

 set(gcf,'Units','centimeters','Position',[0 0 xSize ySize],'PaperUnits','centimeters' ...
     ,'PaperPosition',[0 0 xSize ySize],'PaperSize',[xSize-2.7 ySize-2],'PaperPositionMode','auto')
 
  subplot(4,3,7)
plot(AF0,XX22,'LineWidth',thick,'color','k');
title('Service Workers D','Interpreter','latex','FontSize',10)
axis tight
xlabel('Technology gap, ${A^{F}}/{A^{D}}$','Interpreter','latex','FontSize',10)
set(gca,'FontSize',8,'TickLabelInterpreter','latex')

 set(gcf,'Units','centimeters','Position',[0 0 xSize ySize],'PaperUnits','centimeters' ...
     ,'PaperPosition',[0 0 xSize ySize],'PaperSize',[xSize-2.7 ySize-2],'PaperPositionMode','auto')
 
  subplot(4,3,8)
plot(AF0,XX25,'LineWidth',thick,'color','k');
title('Innovation Difficulty x','Interpreter','latex','FontSize',10)
axis tight
xlabel('Technology gap, ${A^{F}}/{A^{D}}$','Interpreter','latex','FontSize',10)
set(gca,'FontSize',8,'TickLabelInterpreter','latex')

 set(gcf,'Units','centimeters','Position',[0 0 xSize ySize],'PaperUnits','centimeters' ...
     ,'PaperPosition',[0 0 xSize ySize],'PaperSize',[xSize-2.7 ySize-2],'PaperPositionMode','auto')
 
   subplot(4,3,9)
plot(AF0,XX31,'LineWidth',thick,'color','k');
title('Wealth-Income Ratio','Interpreter','latex','FontSize',10)
axis tight
xlabel('Technology gap, ${A^{F}}/{A^{D}}$','Interpreter','latex','FontSize',10)
set(gca,'FontSize',8,'TickLabelInterpreter','latex')

 set(gcf,'Units','centimeters','Position',[0 0 xSize ySize],'PaperUnits','centimeters' ...
     ,'PaperPosition',[0 0 xSize ySize],'PaperSize',[xSize-2.7 ySize-2],'PaperPositionMode','auto')
 
 subplot(4,3,10)
plot(AF0,XX17plot,'LineWidth',thick,'color','k');
title('Skilled/unskilled wage F','Interpreter','latex','FontSize',10)
axis tight
xlabel('Technology gap, ${A^{F}}/{A^{D}}$','Interpreter','latex','FontSize',10)
set(gca,'FontSize',8,'TickLabelInterpreter','latex')

 set(gcf,'Units','centimeters','Position',[0 0 xSize ySize],'PaperUnits','centimeters' ...
     ,'PaperPosition',[0 0 xSize ySize],'PaperSize',[xSize-2.7 ySize-2],'PaperPositionMode','auto')
 
  subplot(4,3,11)
plot(AF0,XX17bplot,'LineWidth',thick,'color','k');
title('Skill Premium F','Interpreter','latex','FontSize',10)
axis tight
xlabel('Technology gap, ${A^{F}}/{A^{D}}$','Interpreter','latex','FontSize',10)
set(gca,'FontSize',8,'TickLabelInterpreter','latex')

 set(gcf,'Units','centimeters','Position',[0 0 xSize ySize],'PaperUnits','centimeters' ...
     ,'PaperPosition',[0 0 xSize ySize],'PaperSize',[xSize-2.7 ySize-2],'PaperPositionMode','auto')
 
   subplot(4,3,12)
plot(AF0,XX19b,'LineWidth',thick,'color','k');
title('Unskilled/Service wage F','Interpreter','latex','FontSize',10)
axis tight
xlabel('Technology gap, ${A^{F}}/{A^{D}}$','Interpreter','latex','FontSize',10)
set(gca,'FontSize',8,'TickLabelInterpreter','latex')

 set(gcf,'Units','centimeters','Position',[0 0 xSize ySize],'PaperUnits','centimeters' ...
     ,'PaperPosition',[0 0 xSize ySize],'PaperSize',[xSize-2.7 ySize-2],'PaperPositionMode','auto')
 
 
print -dpdf comp_stat.pdf
 
